function [Yp,numsv] = hmkp(V,Y)
tic
[n,s] = size(V);

cls = unique(Y);
c = size(cls,1);

sz = zeros(c,1);
idx = cell(c,1);
for i = 1:c
    idx{i}=find(Y==cls(i));
    sz(i)=size(idx{i},1);
end

ix = cell(c,c);
Y = cell(c,c);
P = cell(c,c);

for i = 1:c
    for j = i+1:c
        ix{i,j} = [idx{i};idx{j}];
        Y{i,j} = [ones(sz(i),1);-ones(sz(j),1)];
        P{i,j} = repmat(Y{i,j},1,s).*V(ix{i,j},:);
    end
end

disp('Start optimization...');
tic

ct = 1;
cvx_begin sdp
cvx_solver sedumi
cvx_solver_settings('maxiter',500);
variable t
variable gam(s,nump)
variable M(s,s) symmetric
variable bb(nump,1)
minimize (t);
for i = 1:c
    for j = i+1:c
        [M gam(:,ct);gam(:,ct)' t] >= 0;
        P{i,j}*gam(:,ct)+bb(ct)*Y{i,j}-ones(size(P{i,j},1),1) >= 0;
        ct=ct+1;
    end
end
trace(M) <= n;
%trace(M) - trace(M*TT) <= n;
cvx_end
optime = toc

VM = V*gam+repmat(bb',n,1);
vot = zeros(n,c);

ct = 1;
for i = 1:c
    for j = i+1:c
        vot(:,i) = vot(:,i)+(VM(:,ct)>0);
        vot(:,j) = vot(:,j)+(VM(:,ct)<0);
        ct=ct+1;
    end
end
[~,Yp] = max(vot,[],2);

numsv = sum(del>0);
Yp = Yp-1;
end